Perturbative expansions from Monte Carlo simulations at weak 
coupling: Wilson loops and the static-quark self-energy 
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^ ■ Abstract 
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O 

. Perturbative coefficients for Wilson loops and the static-quark self-energy 

^ : are extracted from Monte Carlo simulations at weak coupling. The lattice 

volumes and couplings are chosen to ensure that the lattice momenta are all 
perturbative. Twisted boundary conditions are used to eliminate the effects of 
lattice zero modes and to suppress nonperturbative finite-volume effects due 
CNl ' to Z{2,) phases. Simulations of the Wilson gluon action are done with both 

periodic and twisted boundary conditions, and over a wide range of lattice 
volumes (from 3^ to 16^) and couplings (from /3 ~ 9 to /3 ~ 60). A high 
precision comparison is made between the simulation data and results from 
finite-volume lattice perturbation theory. The Monte Carlo results are shown 
■ to be in excellent agreement with perturbation theory through second order. 

_ . New results for third-order coefficients for a number of Wilson loops and the 

D , static-quark self-energy are reported. 



I. INTRODUCTION 

Simulations using highly-improved lattice actions have become commonplace in recent 
years. Effective use of these requires perturbative matching calculations for masses, coupling 
constants and currents, among other quantities. Higher-order perturbative calculations for 
these actions are laborious but they are essential in order to obtain precision results for most 
observables. 

An alternative to doing calculations in analytical perturbation theory is to directly mea- 
sure short-distance quantities in Monte Carlo simulations at weak coupling, as proposed in 
Refs. |l|J^. One exploits the fact that the lattice theory on a finite volume enters a per- 
turbative phase at weak coupling. In effect the couplings and volumes in the simulations 
are chosen to ensure that the lattice momenta are all perturbative (up to possible zero 
modes). In this way one can in principle extract perturbative expansions for many quanti- 
ties, by fitting Monte Carlo data for appropriate correlation functions to a power series in 
the coupling. 

This approach has been shown to reproduce analytical results for the first-order mass 
renormalization for Wilson fermions, and the first-order additive self-energy for NRQCD 
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fermions |T|,|^ . In addition, preliminary estimates of some third-order Wilson loop coefficients 
were made in Ref. [H- An extension of this technique to background field calculations was 
considered in Ref. Preliminary work on perturbative simulations of quark actions has 
also been done PJ^. 

The Monte Carlo method as implemented here requires as input from conventional lattice 
perturbation theory only the expansion of the average plaquette, and of the static potential 
(or some other quantity that defines a physical coupling constant), to the desired order, along 
with an estimate of the scale relevant to the quantity of interest. These inputs are necessary 
because we use a renormalization scheme defined through the perturbative expansion of the 
static-quark potential The renormalized coupling ap in this scheme can be extracted 

from measured values of the average plaquette, given its perturbative expansion. 

The method proceeds as follows. Simulations are done at a several different values of j3 
{j3 = 2N/ for SU(iV) gauge theory). At each (3 we use the measured value of the plaquette 
to solve for the value of the renormalized coupling ap{q1^), at the scale g* ^ that is optimal 
for the plaquette . We then run the couplings at each [3 to the scale q* appropriate to the 
quantity of interest, whose expectation values are then fit to a truncated series in ap{q*). 
The fit yields numerical values for the perturbative coefficients. To assess the effects of the 
truncation of the perturbation series at finite order in ap, the fits are done including many 
higher-order terms, beyond the order of interest, but where the fits incorporate constraints 
on the coefficients 0, which are required to lie in a range of values that is consistent with 
a well-behaved perturbative expansion. One can also improve the quality of the results by 
using lower-order coefficients from conventional perturbation theory, if available, in order to 
further constrain the fits to the Monte Carlo data, thereby obtaining more accurate values 
for previously unknown higher-order terms. 

In order to ensure that the lattice momenta are all sufficiently perturbative, simulations 
are done on lattices for which 

fa 

where Aqcd is the QCD scale parameter, a is the lattice spacing, and La is the physical 
length of the lattice. In order to minimize perturbative finite-volume errors we must also 
ensure that the spacing between lattice momenta is small compared to the characteristic 
momentum scale q* associated with the quantity of interest 

In practical simulations the lattices are such that Eq. ([l|) is extremely well satisfied: 
oAqcd ranges from 10~^ to 10~^^ in the analyses in this paper. Moreover the quantities 
studied here all have characteristic scales near the ultraviolet cutoff, hence aq* is indepen- 
dent of a, and Eq. (0) can be easily satisfied. In theories with additional scales, such as quark 
masses m, the characteristic scale q* is proportional to f{ma)/a, where / is some dimension- 
less function, hence aq* is independent of a provided that m is adjusted to keep ma fixed. 
Perturbative finite-volume effects can also be analyzed in detail by running simulations at 
several different volumes. 

In simulations with periodic boundary conditions (PBC) however there are lattice zero 
modes which will violate Eq. (|I]). One must also ensure that nonperturbative finite- volume 
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effects arising from Z{N) phases are sufficiently suppressed. One way to achieve this, using 
a local updating algorithm, is to work on lattices with sufficiently large volumes, where the 
simulation is started with perturbative initial values for the links (e.g., by setting all links 
to the identity). A more powerful approach, analyzed in detail here, is to adopt boundary 
conditions that eliminate zero modes and suppress transitions between different phases. 

It is also desirable to have a more general means for dealing with potential infrared 
problems. The effects of lattice zero modes on most observables is not known and, while 
these could in some cases be accounted for by a numerical fit to the data, it is preferable 
to eliminate these states from the outset. More generally zero modes can significantly alter 
the expected perturbative form of correlation functions. For example, the presence of zero 
modes would pose a significant problem for extracting quark masses, as one could not assume 
the existence of an isolated pole in the perturbative quark propagator. Twisted boundary 
conditions (TBC) |[IOHl2| can be used to eliminate zero modes, and are easily incorporated 
into simulations using existing code for a given action. 

In this work we present a comprehensive study of this Monte Carlo method for extracting 



perturbative quantities [|T3[. We do simulations of SU(3) gauge theory using the Wilson 
gluon action. The evolution of the ap coupling is known for this theory through three- loop 
order which, in principle, allows one to use the method to determine perturbative expansions 
through fourth order. Simulations are done for both PBC and TBC. We show that using 
TBC is an effective means of suppressing nonperturbative finite-volume effects due to Z{3) 
phases, as well as eliminating the effects of lattice zero modes. We also make an extensive 
analysis of perturbative finite-volume effects for TBC. We analyze simulation results for a 
large set of Wilson loops and the static-quark self-energy, over a wide range of couplings 
and lattice volumes. A high-precision comparison is made between the simulation data 
and results from finite-volume lattice perturbation theory. The Monte Carlo results are 
shown to be in excellent agreement with perturbative calculations through second order, 
which are available for these observables for both periodic pi| , p]3| and twisted |]TB[| boundary 
conditions. New results for third-order coefficients for fourteen different Wilson loops and 
the static-quark self-energy are also reported. 

Wilson loops provide a good quantity for a first test of the Monte Carlo method, since 
small loops are relatively insensitive to nonperturbative phases, and have small finite- volume 
errors. The perturbative expansion of small Wilson loops is also relevant to determinations 
of the strong coupling from lattice simulations The calculation of the static-quark self- 
energy is considerably more involved, as it is very sensitive to nonperturbative phases and has 
large perturbative finite-volume corrections. The static-quark self-energy thus represents a 
good prototype for more realistic calculations of other perturbative quantities, such as quark 
masses. The self-energy is also useful for determinations of the 6-quark mass from lattice 
simulations [|17 |. 



The rest of this paper is organized as follows. In Sect. ^ we use PBC to study Wilson 
loops at large /3. Simulations are done on 16^ lattices at nine couplings. We evaluate Wilson 
loops of sizes RxT with R,T < 5. The results are fit to a truncated perturbation series, 
using the renormalized coupling ap evaluated at a scale q^j^ that is appropriate to each 
Wilson loop 0. We explicitly subtract the leading effects of lattice zero modes from the 



Monte Carlo data using an existing analytical calculation [18|. A detailed comparison is 



made with perturbation theory through second order, and estimates are made of the third- 
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order coefficients. 

In Sect, m we review how simulations are done using twisted boundary conditions, and 
their use in analytical perturbation theory. We demonstrate that these boundary conditions 
can be used to virtually eliminate nonperturbative effects due to Z{3) phases even on lattices 
with very small volumes. 

The static-quark self-energy is analyzed in Sect. |^ We extract the self-energy from the 
gauge-invariant Polyakov line, which describes the propagation of a static quark across the 
entire time-extent of the lattice. Other extractions of the self-energy have relied upon large 
Wilson loops | ]I3|p!U[ , and the gauge-fixed quark propagator ||^ (the latter methods could 
prove more efficient in Monte Carlo simulations, as one can limit the propagation time on a 
lattice of a given size). We make a perturbative analysis of the self-energy on finite lattices, 
and show that one can sum leading logarithms in the finite-volume correction. We present 
results of simulations of the self-energy, with runs at nine volumes at each of nine couplings. 
We compare the simulation results with perturbation theory through second order, over the 
whole range of lattice sizes studied here, and make an estimate of the third-order self-energy, 
including an extrapolation to the infinite- volume limit. 

During the course of the calculations described in this paper, we also investigated the 
question of how to design the most efficient calculations through intelligent parameter 
choices, by using the techniques of constrained curve ffiting [Q. In calculations like those in 
this paper, for example, one can choose the couplings (3 in order to minimize the simulation 
cost required to achieve a given precision in the perturbative coefficients. Although the 
Monte Carlo calculations described in this paper were designed before these optimization 
techniques were worked out, we include a description of them here, in an Appendix, for use 
in future calculations. 

Some conclusions and prospects for future work are briefly discussed in Sect. |V|. 



II. WILSON LOOPS 

In this section we analyze simulations of Wilson loops at large (3. The Wilson gauge-field 
action iSwii for SU(3) color is used, where 

Swii[U] = /? E (l - |ReTr U^,{x)) , (3) 

and Ufj,iy is the plaquette. Simulations were done on 16^ lattices at nine couplings. Details 
of the simulation parameters are given in Table |. Periodic boundary conditions were used 
here, in order to make a direct comparison with the ffist- and second-order perturbative 
coefficients calculated on finite lattices by Heller and Karsch |[T^. We will also verify that 
the simulations reproduce the effects of lattice zero modes to leading order. 
We analyze the logarithm of the Wilson loop 

~ 2{r\t) ^""'^ = T.^n apilkr)^ (4) 

using a renormalized coupling ap that is determined from measured values of the plaquette, 
according to MM 
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TABLE I. Simulation parameters for Wilson loop measurements. These simulations were all 
done on 16^ lattices with periodic boundary conditions. The lattice coupling (3 for each simulation 
was determined from the bare coupling aiat by /3 = 6/ (47raiat ) • The measured values of the 
average plaquette are shown, along with the renormalized couplings ap(3.40/a) and scale masses 
Ap extracted from Eqs. (^) and (|^). Ten configurations were skipped between measurements in all 
cases, and the observables were computed by binning the measurements in bin sizes of 50, which 
resulted in negligible autocorrelations at all couplings. 



aiat 




Measurements 


(iRe-TrC/n) 


ap(3.40/a) 


aAp 


0.010 


47.746 


2320 


0.957542(1) 


0.01049 


5.47 X 10-23 


0.015 


31.831 


2459 


0.935857(2) 


0.01614 


8.63 X 10-^^ 


0.020 


23.873 


2112 


0.913829(3) 


0.02209 


1.05 X 10-1° 


0.025 


19.099 


1558 


0.891441(3) 


0.02839 


2.94 X 10"^ 


0.030 


15.915 


860 


0.868599(9) 


0.03510 


1.25 X 10-^ 


0.035 


13.642 


746 


0.845305(8) 


0.04225 


1.81 X 10^5 


0.040 


11.937 


748 


0.821472(9) 


0.04992 


1.34 X 10-^ 


0.045 


10.610 


500 


0.797038(11) 


0.05819 


6.40 X lO"'' 


0.050 


9.459 


500 


0.771872(11) 


0.06719 


2.24 X IQ-^ 



-In 1^1 



47r 



ap(3.40/a) [1 - 1.1909 ap] . 



(5) 



The coupling ap is defined such that the logarithm of the plaquette has no third- or higher- 
order terms in its perturbative expansion Other quantities, of course, do have higher- 
order terms when expressed as a series in ap. One can express ap(3.40/a) as a series in the 
bare lattice coupling aiat, using the third-order expansion of the plaquette given in Ref. Il2( 



ap(3.40/a) 



aiat + 4.564 afat + 28.566 afat + 0(a^ 



lat.* 



(6) 



The large coefficients in this expansion are an artifact of aiat; using ap eliminates large 
renormalizations of the bare coupling. We also note that the logarithm of the Wilson loop is 
better behaved perturbatively than the Wilson loop itself, due to the exponentiation of the 
perturbative perimeter law; this is also why we have divided by the perimeter in defining 
the perturbative coefficients in Eq. (^). 

The perturbation series for each Wilson loop is evaluated using the renormalized coupling 
at a scale g|j determined according to the procedure of Ref. 0; the scale corresponds to the 
typical momentum carried by a gluon in the leading-order diagram for a given quantity. The 
scales are given in Table ||. The couplings at these scales are evaluated by first measuring 
the average plaquette in the simulation, and solving Eq. (|]) for ap(3.40/a). We then evolve 
to the scale appropriate to the quantity under consideration using the universal second-order 
beta function, plus the third-order term for ap, with 



ap{q) 



An 



PI 



X 



In 



1\2 



/3? 



(7) 
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TABLE IL Scale parameters for the couplings for various Wilson loops. 



Loop '^Qrt Loop (^Q*RT 

1 X 2 3i)7 2 x 5 2.46 

1 X 3 3.01 3 x 3 2.45 

1 X 4 2.96 3 x 4 2.38 

1 X 5 2.95 3 x 5 2.35 

2 x 2 2.65 4 x 4 2.30 
2 x 3 2.56 4 x 5 2.27 
2 x 4 2.49 5 x 5 2.23 



For our quenched simulations (3o = 11, (3i = 102, and (32,p = /52MS + where (32MS ~ 
2857/2 is the third beta function coefficient in the MS scheme, and where Bp = —147.57 
can be obtained from existing three- loop calculations PD|P^ , as described in Ref. The 
values of ap(3.40/a) and Ap for our simulations are given in Table |. 

One can convert a perturbative expansion in a\s,t to one in ap{q), at the scale appropriate 
to a particular quantity, through third order, using 



aiat = Oip{q) - aj,{q) 
+ ap(g) 



A-K \aq I 



+ 4.702 




+ 4.702 




7.841 ) +0{a 



(8) 



This connection can be obtained from the third-order expansion of the plaquette in the bare 
coupling which can be used to solve for aiat in terms of ap(3.40/a), given its definition 



which can be used to solve for ai^x in terms of ap(3.40/a) 
one can then use a perturbative expansion of the evolution equations 



to 

extends the second-order connection 



in Eq. (| 

eliminate ttp(3.40/a) in favor of ap{q). Equation 
between aiat and Oip{q) given in Ref. |0. 

The ffist- and second-order perturbative coefficients for the Wilson loop were computed 
by Heller and Karsch for an expansion in the bare lattice coupling. We convert this 
expansion to a series in Cip{q*pj) using Eq. (H). These perturbative calculations were done 
on finite lattices with periodic boundary conditions, neglecting the contribution of lattice 
zero modes. Our simulations of the Wilson loops were also done with PBC but do contain 
the effects of zero modes. In Sect. |IT1| we will consider simulations using twisted boundary 
conditions to eliminate these states. For Wilson loops however we can make use of an 
analytical calculation of the zero mode piece c^^''° of the ffist-order coefficient in Eq. @, due 



to Coste et al. [jT8 



47r(i?T)2 
9{R + T)V'' 



(9) 



where V is the lattice volume. We will use this expression to explicitly subtract the leading- 
order effects of the zero modes from the Monte Carlo data. In the following we will show 
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the Monte Carlo data after first making this zero mode subtraction, unless explicitly noted 
otherwise; hereafter we will use c„ to denote the (finite-volume) coefficients without zero 
mode contributions. 

We present Monte Carlo results for the 5x5 Wilson loop in Fig. |I|, in terms of the 
quantity 



MC 



ap{q%. 



2{R + T) 



R,T 



(10) 



which should exhibit the limit ki — >■ Ci as ap — >■ 0. We extract estimates of the perturbative 
coefficients from the Monte Carlo data by fitting the results to the series expansion Eq. 
(I) where, to begin with, we treat all coefficients c„>i as unknown. We will compare the fit 
values for ci and C2 with the results from analytical perturbation theory, which provides a 
stringent test of the Monte Carlo method. 

An important aspect of the fitting procedure is how to reliably account for the systematic 
error arising from the truncation of the fit function at a finite order in ap. Including too 
few terms in the expansion in ap results in a poor fit to the data, while including too many 
higher-order terms results in very poorly constrained values for the lowest-order coefficients 
which should, in fact, make the dominant contributions to the data. This situation can 
be remedied by incorporating constraints on the coefficients, which are required to lie in a 
range of values that is consistent with our expectation that the perturbative expansion is 
well behaved. We do this in practice by using conventional least-squares fitting routines, 
where the is augmented according to: 



^ ^2 



X^(Cn) - xL,{Cn) = X\C^) + E '-^^^^-P^, (11) 

which tends to constrain the fit values for the c„ to the interval c„ ± o"„ . This approach can 
be motivated by Bayesian statistical analysis [§. 

If perturbation theory is reliable we expect the coefficients c„ to be of 0(1). We performed 
least-squares fits to Eq. (|), minimizing x^^g with c„ = and o"„ = 5 for the first five orders 
in the expansion. The dashed line in Fig. |I] shows the results of the fit for the 5x5 Wilson 
loop; note that the curvature in the Monte Carlo data shows the sensitivity of the 
simulations to the third-order term in the perturbative series. The quality of the fits is very 
good, with Q-values in excess of 50%. 

The measured values of Ci and C2 are in excellent agreement with perturbation theory, 
as shown in Table PTI] , with an accuracy of a few parts in 10^ for the first-order coefficients 
and a few parts in 10^ for the second-order coefficients. The third-order coefficient can also 
be resolved, here with almost no input from analytical perturbation theory. The fit values 
are very stable to changes in the values of c„ and ct„ used in Eq. (pjj). 

Note that if the Monte Carlo data are fit with ci constrained to its value from pertur- 
bation theory, then the errors on C2 are reduced by a factor of about three, with fit values 
in agreement with perturbation theory within the reduced errors. Similarly, we obtain more 
accurate results for C3 by fitting the Monte Carlo data with ci and C2 constrained to their 
perturbative values. We did fits to Eq. (^) for the next three orders in the perturbative 
expansion, minimizing x^^g using c„ = and a„ = 5 for = 3, 4, 5. The results for C3 are 



given in Table IV, where the fit errors are seen to be about 10%. 
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TABLE III. Monte Carlo results for the first three perturbative coefficients for selected Wilson 
loops {0^23)- The results were obtained from a simultaneous fit to the coefficients, as discussed in 
the text. The first- and second-order coefficients from perturbation theory for the same size lattice 
are also shown (0^2^). The effects of zero modes are not included in the perturbation theory values, 
and were removed at leading order from the simulation data. Note that more accurate results for 

for the full set of Wilson loops are ] 
constrained to their perturbative values. 



for the full set of Wilson loops are given in Table IV, where the fits are done with ci and C2 



Loop 


^1 


^1 


„MC 
^2 


'^2 




1 X 2 


1.2037(2) 


1.2039 


-1.244(16) 


-1.260 


0.0(5) 


1 X 3 


1.2587(2) 


1.2589 


-1.185(19) 


-1.198 


0.4(5) 


2x2 


1.4337(2) 


1.4338 


-1.312(19) 


-1.323 


1.1(5) 


3x3 


1.6089(3) 


1.6089 


-1.218(24) 


-1.217 


2.5(6) 


4x4 


1.7067(4) 


1.7067 


-1.213(29) 


-1.210 


3.4(6) 


5x5 


1.7693(6) 


1.7690 


-1.201(40) 


-1.177 


4.3(7) 



TABLE IV. Monte Carlo results for C3, where the first- and second-order coefficients are con- 
strained to their values from perturbation theory. 



Loop 



Loop 



1 X 2 
1 X 3 
1 X 4 
1 X 5 
2x2 
2x3 
2x4 



0.43(9) 

0.66(11) 

0.84(12) 

0.94(14) 

1.41(11) 

1.91(13) 

2.28(15) 



2x5 
3x3 
3x4 
3x5 
4x4 
4x5 
5x5 



2.52(17) 
2.53(15) 
2.98(17) 
3.26(19) 
3.40(19) 
3.71(21) 
3.91(23) 



It is also interesting to verify that the simulations reproduce the leading effects of the 
zero modes. A convenient way to visualize these effects is to plot the quantity 



— ^2 



aUq_ 



R,T) 



2{R + T) 



(12) 



where the first-order coefficient is set to its perturbative value. We plot K2 for the 5x5 
loop in Fig. ^, after subtracting the leading-order zero mode term from the Monte Carlo 
data. We see that the data reproduce the second-order coefficient from perturbation theory, 
with ^2 — > C2 as ap — > 0. In Fig. ^ we plot K2 but where the leading-order zero mode is 
not subtracted from the data. We see evidence of singular behavior in K2 at small coupling, 
indicating that the first-order term is not completely removed from the Monte Carlo data 
when the zero mode component is not treated. The dashed line in Fig. ^ shows the results 
of a fit to Eq. (^), taking account of the leading zero mode contribution, where the term 
c^i^°/ap is included in the fit line. This shows explicitly that the Monte Carlo data are 
sufficiently accurate to reveal the small contribution from the zero modes, at sufficiently 
small couplings. 

We also present results for the residual 
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^3 



MC _ 

~ a%{q*R,T) 



(13) 



for the 5x5 loop in Fig. ^, which is convenient for visuahzing the sensitivity of the Monte 
Carlo data to the third- and fourth-order terms. The statistical errors in the Monte Carlo 
data are too large to resolve C4, although the best fits suggest that C4 is of the same order 
as C3 for all the Wilson loops analyzed here. 

A potential complication in our analysis of C3 is that we have only corrected the Wilson 
loop data for the effects of zero modes to first order. However we expect that the leading 
contribution from zero modes that remains is of Ola^lRT)"^ /V) which, given the lattice 
volume and the range of couplings analyzed here, should only be comparable to terms of 
0(ap). In fact there is no visible effect of zero modes beyond first order, within statistical 
errors; this would should up as singular behavior in ^3 at small ap (compare Fig. |^ for 
with the two plots of K2 in Figs. || and |]). We also note that while we have extracted values 
for C3 on a finite lattice, the volume is large enough that the results should give a good 
approximation to the coefficients on an infinite lattice (with the corrections expected to be 
oiOil/V)). 

A determination of higher order terms in the expansion of the 1x1 and 2x2 Wilson loops 
has also been made in Ref . , using numerical simulations of the Langevin equations, where 



a perturbative expansion in the bare lattice coupling is applied to the evolution equations 
themselves. The results are presented in Ref. as an expansion of the Wilson loops in 
powers of the bare lattice coupling, Wr^t = 1 — SnCnC^M- We can convert our third-order 
result for the expansion of IniWR^T) in ap{q}^j.), to an expansion of Wr^t in aiat, by using 
the inverse of Eq. (H). We find C3 = 0.3 ± 0.9 for the 2x2 Wilson loop, in agreement with 



the result C3 = 0.0 ± 0.9 reported in Ref. |24 



We note parenthetically that the expansion of the Wilson loop itself is very poorly 
convergent, and that the vanishingly small value of C3 for the 2x2 loop is accidental. One 
finds very large expansion coefficients for other Wilson loops. These expansions are tamed 
by taking the logarithm, and expressing the series in ap{q*). For example, our results give 
C3/C1 = 76.1 ± 0.1 for the 5x5 loop, with this large value arising almost entirely from 
the exponentiation of the perturbative perimeter law, and the renormalization of the bare 
coupling (compare with c^/ci = 2.2±0.1 for the logarithm of the 5x5 loop). In Sect. [V we 
show that the third-order expansion of the static-quark self-energy is also very reasonable 
when expressed in terms of ap{q*), but is very poorly convergent when expressed in terms 
of aiat- 



III. TWISTED BOUNDARY CONDITIONS 
A. Formalism 

The analysis of the preceding section shows that simulations at large /3 can be used to 
make accurate determinations of perturbative quantities at higher orders than have been 
achieved using conventional perturbation theory. However, as discussed in Sect. |, simu- 
lations with periodic boundary conditions (PBC) are subject to the effects of lattice zero 
modes, and a more convenient method for dealing with potential infrared problems in more 
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general situations is required. Fortunately lattice zero modes can be completely eliminated 
by using twisted boundary conditions (TBC). We also find that TBC significantly reduce 
nonperturbative finite- volume effects due to Z{N) phases. 

Twisted boundary conditions [p!0H l^ for the link fields resemble a gauge transformation 
on fields which cross through selected lattice boundaries 

u^{x + LO) = n,Ua{x)nl (14) 

where we take the twist matrices Qi, to be constant. A twist must be applied to at least 
two boundaries fi and u, as the transformation matrix Qi, for a single boundary can be 
completely eliminated by a field redefinition. The requirement that Ua{x + Lfi + Lz>), which 
can be connected to Ua{x) by crossing the two lattice boundaries in different orders, be 
single-valued implies that the twist matrices must satisfy the algebra 

n^n, = r]n,n^, vez^N), (15) 

where we consider the general SU(iV) gauge theory in most of this section. A pair of twist 
matrices generates a multiplication table that forms a discrete subgroup of SU(A^). In 
particular |ll|] 

^A. ^ (_1)A^-1/, (16) 

where / is the unit matrix. 

The Wilson gluon action for TBC is written in terms of link variables in the usual way, 
Eq. (^. The matrices A(a;) that generate gauge transformations 

U^ix) A{x)U^ix)A\x + fi) (17) 

are subject to the same TBC as the links (Eq. (0)). One consequence of this is that the 
Polyakov line in a twisted direction ^ must have an additional factor of the corresponding 
twist matrix, if it is to be gauge- invariant: 

= {u^{x)u^{x + fi)... u^{x + LjX) X n;) . (is) 

The only zero-action fields are pure-gauge configurations U^{x) = A(x)A"l'(x + fi) [0, in- 
cluding possible Z{N) phases. The action also possesses a discrete symmetry 

u^{x) ^iU^nl (19) 



where fli = fi^, Vl^, etc. Note that Eq. ([T9|) is not a gauge transformation, because 

A{x) = Q does does not satisfy TBC ||TT|. This symmetry implies that the Polyakov lines 
in the twisted directions have zero expectation value, even in the perturbative phase of the 
theory. 

We have done SU(3) simulations with twisted boundary conditions across two "spatial" 
boundaries x and y 

^x^y = V^y^x (J^xy), (20) 
where r] = e^'^*/^, and across all three spatial boundaries x, y and z 
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(21) 



We refer to these two CcLSGS cLS Txy and Txyz boundary conditions, respectively. 

Explicit representations of the twist matrices are not needed since, as shown in Ref. 
Ill , one can absorb them by a field redefinition of the link variables, leaving only phase 



factors 7] and rj* multiplying the plaquettes at the corners of the twisted planes. We have 
done simulations with this field redefinition, and also using an explicit representation of the 
boundary conditions Eq. (0) using the matrices (for SU(3)) 



(22) 





"0 


1 


0" 

















1 




1 







1 











g27rj/3 



and 





1 

e^^^/^ 



-2ni/3 






(23) 



The boundary conditions Eq. (|T^ lead to an unusual quantization of the lattice momen- 
tum modes, as well as removing the zero modes. Making the usual substitution 



JgaAf^(x) 



(24) 



the boundary conditions take the form A^(x + LO) = Q^Af^{x)Ql. A twisted plane wave 
basis is used to Fourier analyze the fields 



^.(^) = ^Ex.r,e^'-(^+^^^4(fc), 



(25) 



where, in order to obey the boundary conditions, the matrices must satisfy the algebra 

n^YkQl = e^'^-^Ffc, (26) 

and where Xk enforces a constraint on the mode sum, to be developed below (see Eqs. ( pO[ ) 
and (|31|)). 

The quantization conditions follow by iterating Eq. (^) times and using Eq. (jTB]). 
One finds that momenta in twisted directions are quantized as if the SU(A^) fields live on a 
lattice of length Lx N, rather than the actual length L (although some modes are excluded) 



2tt 

1 

LN 

271 



L 



twisted direction, 



periodic direction. 



(27) 



The extra momentum degrees-of-freedom come about because the color structure of each 



mode is unique, up to a phase. Substituting F^ 
convenient choice of phase [pT| , 



into Eq. 



one finds, with a 
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Tk = n-^'y n"^- ^2 (28) 
These matrices are orthonormal under the trace 



-Tr (rl,r,) = 1 1 '[<y = (29) 

\ ^ "J \0 otherwise. ^ ^ 

Since the fields must be traceless one finds that a set of modes, including the zero modes, 
are excluded 

_ f i{n^ = ny = (mod A^), . . 

\ 1 otherwise. 

In the case of Txyz boundary conditions, a further constraint emerges from Eq. 
Hence there is a factor of A^^ — 1 more momentum modes with TBC, each of which has a 



single color degree-of-freedom, which is exactly the number of independent colors that one 
has for each momentum mode with PBC 



At tree-level the twisted gluon propagator in momentum space has the structure JIT 



1 , 



{A^{k)A,ik'))g=o = lVNxkV'~^^'''''^6,,k>D^,{k), (32) 
where 

{k', k) = n'^n^ + riyUy + {n^ + ny){n^ + n'y). (33) 
For the Wilson action in Lorentz gauges one has 



k k 
k^ 



(34) 



with k^j, = 2sin(|A;^) and k"^ = J2x ^, 



B. Suppression of Z{N) phases 

The Polyakov line along an untwisted direction is an order parameter for the Z{N) 
degenerate vacua of the lattice theory, which correspond to the invariance of the Wilson 
action under the transformation 

f/^(x) — » riU^{x), Wx 3 X ■ fi = constant. (35) 

The Polyakov line is also sensitive to the formation of domains between different Z{N) 
phases. These nonperturbative effects must be suppressed if one is to use Monte Carlo 
simulations at large (3 to extract perturbative quantities. In particular we will use simulation 
results for the Polyakov line itself to obtain the perturbative self-energy of a static quark. 
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In order to suppress nonperturbative Z{N) phases we start the simulation with all links 
initialized to a "cold start," t/^ = /. The probability of making a transition to an nontrivial 
Z{N) phase in a local updating algorithm can then be reduced by working at sufficiently 
large lattice volumes. In fact our results for Wilson loops on 16^ lattices with PBC, presented 
in Sect. |T|, are in excellent agreement with finite- volume perturbation theory. However we 
find that nonperturbative Z{N) phases are generated frequently on small lattices when PBC 
are used, and this occurs even at extremely large (3. On the other hand, we find that using 
TEC leads to a dramatic suppression of these effects compared to PBC, on lattices of the 
same size. 

We illustrate the effects of Z{3) phases with simulation results for 4^ lattices at /3 = 9. 
We show run time histories and scatter plots of the real and imaginary parts of the Polyakov 
line along an untwisted direction (hereafter taken to be the "temporal" direction t), where 

Pt{L) = ^ E ReTr ^ fl f/,(x)^ , (36) 

for a lattice of volume L^. Results for PBC are shown in Fig. ^, for Txy boundary conditions 
in Fig. ^ and for Txyz boundary conditions in Fig. ^ We see that nonperturbative Z{3) 
phases and domains render simulations with PBC useless for extracting perturbative quanti- 
ties on small lattices. We also see that twisted boundary conditions create a barrier between 
Z(3) phases, and that transitions between these phases and are essentially eliminated with 
Txyz boundary conditions (with no tunneling events observed in millions of updates in the 



range of (3 values considered here). In Sect. |^ we show that the remaining finite- volume 



effects on lattices with Txyz boundary conditions are very well described by perturbation 
theory for /3 ^ 9, even on volumes as small as 3^. 



IV. STATIC-QUARK SELF-ENERGY 

A. Perturbation theory 

In this section we consider the perturbative expansion of the self-energy Eq of a static 
quark. We extract the self-energy from the gauge-invariant Polyakov line Pt along an un- 
twisted direction, which describes the propagation of a static quark across the entire time- 



extent of the lattice. One could also obtain the self-energy from large Wilson loops [|T5|,|T9 
or from the gauge-fixed static-quark propagator This study however represents a good 
prototype for calculations of other more realistic perturbative quantities, such as quark 
masses. 

We ffist define the self-energy Eq{L) on a finite lattice according to 

aE^{L) = ~\n{Pt{L)) (37) 

where, for comparison with our simulation results in the next subsection, we consider lattices 
with equal lengths L along all sides. One then obtains the infinite-volume self-energy Eq by 
taking the limit, 

= Eo(L ^ oo). (38) 
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We analyze the tadpole-improved self-energy. This is obtained by dividing the links 
in the Polyakov line by a mean field Mq, Ufj_{x) U^{x)/uq. Hence the tadpole-improved 
self-energy is related to the unimproved self-energy by the addition of In(Mo). We use the 
average plaquette to estimate the mean-field: 

uo = {UuY'\ (39) 

The expansion of the self-energy to second order was computed in perturbation theory 
according to Eqs. ( p7|) and (|38D by Heller and Karsch, for an expansion in the bare coupling, 
with the result [0,|1|] 

^^unimp ^ 2.1173aiat + 11.152< + 0«) (40) 

for the unimproved self-energy. Hereafter we consider only the tadpole-improved self-energy, 
which we denote by Eq. We convert Eq. (50) to an expansion in the renormalized coupling 



at the appropriate scale using Eq. (||): 

aEo = 1.0701 ap(g^J + 0.117a| + 0(a|), q^^ = 0.84/a, (41) 

where Eq. (|]), with couplings evolved to q^^, provides the tadpole subtraction. 

In the next subsection we will compare Monte Carlo data for the self-energy with results 
from analytical perturbation theory on finite lattices. In order to extract the infinite-volume 
self-energy Eq from the Monte Carlo simulations we must also make an extrapolation of 
measurements of Eq{L) done on finite volumes. We can gain some insight into the nature 
of the perturbative finite-volume corrections from some analytical considerations. 

We define perturbative coefficients c„(L) on a finite lattice according to 

aEo(L) = ^c„(L)a^(g|;J. (42) 

n 

For TBC the first-order term is given by 

c,{L) = ^Y.XkD,,{k, = Q,k)-^ [TBC], (43) 

k 

to be compared with a calculation for PBC, ignoring the contribution from zero modes 

ci{L)= ^ E ^44(fc4 = 0, fc) - - [PBC]. (44) 

The constant 7r/3 that is subtracted from the momentum sums in the above expressions 
is the value of the one-loop tadpole- improvement counterterm, neglecting its very weak 
dependence on the lattice volume. We remind the reader that the mode sums in Eqs. ( ^31) 
and ( p4|) are different, due to the different quantization of the momentum components along 
the twisted and periodic directions. The two sets of boundary conditions yield identical 
results in the infinite- volume limit, where the color factor A^^ — 1 emerges in the case of 
TBC because Xk averages to A^^ — 1 over infinitesimal momentum intervals. 

Results for Ci{L) for the three boundary conditions are presented in Fig. H, which shows 
that finite- volume effects are reduced with TBC, as suggested by Eq. (|27|). As is evident 
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from the plots, finite-volume corrections are very well parameterized by a simple linear form 
in 1/L, 

ci(L)=ci-Xi^ + 0(-^), (45) 



where Ci = Ci{L = oo). One can evaluate Xi numerically from Eqs. (|43| ) and (^^ with the 
results Xi ^ 1.891 (PBC), 0.254 (Txy), and 0.771 (Txyz). 

We expect the finite- volume correction to run with a coupling ap{q1), evaluated at an 
infrared scale q1 that is set by the box size: 



aEo{L) = aEo-X,^ + O r^,^ , gl oc f (46) 




A physical interpretation of this functional form for Ci(L) is that the static quark experiences 
a perturbative Coulomb interaction with its images in the walls of the lattice. The different 
values of the coefficient Xi for different boundary conditions also has a natural interpretation 
in this picture: TBC reduce finite-volume effects by effectively putting the image charges 
further away from the source charge. 

Having established Eq. (|^), one can deduce logarithms in L in the self-energy at higher 
orders. For example, at second order one has 

C2{L) =C2-^{X, + Y, ln{L')) + O (^-^^ , (47) 

where C2 = C2(-C/ = oo), and 

y. = x,|. (48) 



This follows from an expansion of the running coupling in Eq. (46), to second order in the 
coupling at a reference scale, such as ap{qEo)- One can explicitly isolate the logarithm in 
the second order coefficient using existing perturbative calculations, which were done long 
ago by Heller and Karsch in the case of PBC |]T^ , and which have also recently been done in 
Ref. |jl6|] for TBC. Results of the perturbative calculations for the three boundary conditions 



are plotted in Fig. y over a range of lattice sizes. The dashed lines in Fig. y show fits to Eq. 
(^), where C2 is constrained to the correct value; the fits are in excellent agreement with Eq. 
(^). Note that the curvature in the results for C2{L) reveals the presence of the logarithm, 
particulary in the case of PBC and Txyz boundary conditions, where the logarithm makes 
a significant contribution at the lattice sizes shown in Fig. |^. 

In the next subsection we will use Eq. (^) to deduce the form of the logarithms in L in 
the third order self-energy, which will help to constrain the infinite-volume extrapolation of 
the Monte Carlo data. We note that one should similarly be able to determine the leading 
logarithms in the finite-volume corrections to other quantities, which should likewise prove 
useful in Monte Carlo determinations of their perturbative expansions. 
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TABLE V. Simulation parameters for static self-energy measurements with Txyz twisted 
boundary conditions. At each /3 simulations were run on nine volumes, L = [3, 11] inclusive. 
The measured average plaquette on the lattices with L = 10 are given, along with the scale mass 
Ap computed from Eqs. (^) and (^. The bare coupling is shown along the renormalized coupling 
ap{q%) evaluated at the scale appropriate to the self-energy with tadpole renormalization. 



p 


(iRe-TrC/n) 


aAp 


"lat 


ap(0.84/a) 


60.0 


0.966311(1) 


2.57 X 10-29 


0.008 


0.00843 


23.8 


0.913831(4) 


1.05 X 10"^° 


0.020 


0.02338 


19.0 


0.891415(5) 


2.95 X 10"*^ 


0.025 


0.03058 


16.0 


0.869332(6) 


1.13 X 10-^ 


0.030 


0.03824 


13.6 


0.845310(7) 


1.81 X 10-5 


0.035 


0.04731 


12.0 


0.822493(8) 


1.25 X 10"^ 


0.040 


0.05676 


10.6 


0.797025(10) 


6.41 X 10"^ 


0.045 


0.06844 


9.5 


0.771866(11) 


2.24 X 10-3 


0.050 


0.08138 


9.0 


0.756142(13) 


4.29 X 10-3 


0.053 


0.09032 



B. Self-energy from Monte Carlo simulations 

We measured the static-quark self-energy in simulations done with Txyz twisted bound- 
ary conditions at nine couplings. The simulation parameters are given in Table [V]. Simula- 
tions were run on nine volumes L^, L = [3, 11] inclusive, at each of the nine couplings, for 
a total of 81 lattices. The number of measurements made on each volume, at all couplings 
except (3 = 60, were as follows: 2000 measurements for the lattices with L = [3, 6] inclu- 
sive, 1500 measurements for L = 7, 1200 for L = 8, 800 for L = 9, 600 for L = 10, and 
400 for L = 11 (ten times as many measurements were made on each volume at /3 = 60). 
One hundred configurations were skipped between measurements at all couplings, except 
at /3 = 60, where ten configurations were skipped between measurements. The observables 
were computed by binning the measurements in bin sizes of one hundred, which resulted in 
negligible autocorrelations at all couplings. The static energy and its error were computed 
from the binned ensembles using a standard jackknife analysis. 

To demonstrate the reliability of the Monte Carlo method, we first use the simulation 



results to estimate the first- and second-order perturbative coefficients. In Fig. [1^ we plot 
the quantity 

KnL) = aE^^{L)/ap{q%J (49) 

versus ap{qEo), for all values of L. The dashed fines show the results of least-squares fits 
to Eq. (1^), minimizing xing (PI) ) using c„ = and a„ = 5 for the first five orders in 
the expansion. The quality of the fits is very good in most cases, with Q-values typically in 
excess of about 20%, although the lowest Q- value in the fits is 3%. 

We show the Monte Carlo results for Ci(L) in Fig. |TT|, where they are compared with 
finite- volume perturbation theory for Txyz boundary conditions, Eq. (|4^). The data agree 
with perturbation theory within errors of only a few parts in 10^. 
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Monte Carlo results for C2{L) are shown in Fig. |T^. In this case the fits to Eq. (^21) 
were done with Ci(L) constrained to its perturbative value. We see that the Monte Carlo 
simulations also reproduce the results of second order Txyz perturbation theory over 
the full range of lattice sizes, within errors that are as small as a few parts in 10^ at several 
volumes. 

The third-order term in the self-energy is not known from conventional perturbation 
theory. We use our simulation results to estimate c^^L), by redoing fits to Eq. with both 
Ci(L) and C2{L) constrained to their perturbative values. In order to determine the value of 
C3 at infinite volume, we must also account for the systematic error due to the extrapolation 
from finite lattices. The leading finite-volume corrections come from logarithms in L, which 
can be determined using renormalization-group arguments (cf. Eqs. (^6D -(^8l)). Making 
use of these constraints considerably improves the accuracy of the extrapolation to infinite 
volume. 

We show the Monte Carlo data for the third-order coefficient as a function of lattice size 
in Fig. after subtracting the logarithms at 0{1/L); the data are presented in terms of 
the residual 

6c,{L) ^ c,{L) + 1 In^ (^] + 1 In ( ^] . (50) 



The scale length Lq in the leading logarithm is determined from second-order perturbation 
theory where, according to Eq. (^Tf), Lq = exp(— X2/212)- We evaluate the scale length from 
a fit to the Txyz perturbation theory results illustrated in Fig. |^, which gives Lq ^ 0.45 for 
the expansion in ap{q%^). 

The Monte Carlo data are consistent with the logarithms in Eq. (0), which are found 
to dominate the extrapolation to the infinite- volume limit, even from these relatively small 
lattices. To extract the infinite-volume coefficient C3 we fit the remaining finite-volume 
corrections to the form 

1 1 ^ /L2\ 

6Cs{L) = C3 + Pl,o- +T.T7^I1 Pm,n 72 • (51) 

^ m>2^ n=0 \^^J 

Figure shows the results of a fit to Eq. (^) where xLg minimized using C3 = p^.n = 0, 
and = „ = 4, for m = 2, 3 and = 0, 1, 2 (and for pi^). The fit yields 

c^c ^ 3 5g ^ 0.50 (infinite- volume limit). (52) 

Changing the order of the expansion in 1/L in Eq. (|5lD makes litle change in the fit value for 
C3. We conclude from these results that renormalized perturbation theory for the tadpole- 
improved self-energy is well behaved through third order, with the data in Fig. ^ clearly 
demonstrating that C3 is of 0(1). 

An estimate of the third-order term in the expansion of the unimproved self-energy, in 
the bare lattice coupling (Eq. (^0]) ), has recently been reported using numerical simulations 
of the Langevin equations |T9|. We can convert our result for C3 from an expansion in 
a;p(0.84/a) to one in aiat, using the inverse of Eq. (||). We find (^\ix. = 86.6 ± 0.5 (without 
tadpole improvement), in agreement with the value C3jat = 86.2 ± 0.6 reported in Ref. 
||T9| . The bare coupling is clearly a very poor expansion parameter, with 96% of C3 iat being 
absorbed by renormalization when a physical coupling is used. 
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We note that the simulations in Ref. were done on much larger lattices than were 
used here. We were able to extract C3 from smaller lattices because the leading finite- volume 
corrections were identified using renormalization-group methods. This determination of the 
third-order self-energy involved a modest computational effort. The entire set of simulations 
in the present analysis required only the equivalent of about 150 days of running on a single 
1 GHz processor. 



V. SUMMARY AND OUTLOOK 

The results presented here demonstrate that higher-order perturbative expansions are 
accessible in Monte Carlo simulations at large (3. An extensive theoretical analysis was 
presented together with the results of numerical simulations of a large set of Wilson loops 
and the static-quark self-energy. Twisted boundary conditions were investigated as a means 
of eliminating zero modes and suppressing nonperturbative finite-volume artifacts, and an 
extensive analysis of perturbative finite-volume corrections was made. Wilson loops pro- 
vided a good quantity for a first test of the Monte Carlo method, since small loops are 
relatively insensitive to finite-volume effects. The calculation of the static-quark self-energy 
was considerably more involved, as it is very sensitive to nonperturbative phases and has 
large perturbative finite-volume corrections. The static-quark self-energy thus represents a 
good prototype for calculations of other more realistic perturbative quantities, such as quark 
masses. 

The simulation results were shown to reproduce perturbation theory on finite lattices 
through second order to high precision, over a wide range of lattice sizes and couplings. 
Monte Carlo results for the fourteen smallest Wilson loops were found to agree with per- 
turbation theory within the errors, with an accuracy of a few parts in lO'^ for the first-order 
coefficients, and a few parts in 10^ for the second-order terms. The Monte Carlo results for 
the static-quark self-energy were found to agree with finite- volume perturbation theory over 
the full range of lattice sizes analyzed here, with an accuracy of a few parts in 10^ at first 
order, and a few parts in 10^ at second order. This precision was achieved with relatively 
little computational effort. New estimates of third-order terms for the Wilson loops and 
the static-quark self-energy were obtained to about 10% accuracy. Renormalization-group 
arguments were used to improve the quality of the extrapolation of the self-energy to infinite 
volume. The results demonstrate that renormalized perturbation theory for Wilson loops 
and the self-energy is well behaved through third order. 

These methods can be directly applied to improved gluon actions, and can be extended 
to quark actions. We have done some work on large j3 simulations for fermions in the 
nonrelativistic formulation of QCD, extending the preliminary studies reported in Refs. [|l],3. 
We find that simulations of the additive energy and multiplicative mass renormalization 
reproduce results of one-loop perturbation theory, and can resolve the second-order terms 
in the expansion of these quantities, over a wide range of bare quark masses Further 
work in this direction is in progress. 
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APPENDIX: 

DESIGNING OPTIMIZED MONTE CARLO SIMULATIONS 

In this Appendix we consider the question of how to design the most efficient calculations 
through intelligent parameter choices, by using the techniques of constrained curve fitting 
0. As discussed in the Introduction, one objective of this analysis, in the context of short- 
distance Monte Carlo simulations like those in this paper, is to choose the couplings a 
for the simulations so as to minimize the cost required to achieve a given precision in the 
perturbative coefficients. 

As described in Ref. ^ and in Sect. fT|, the effects of truncation errors in fitting power 
series to Monte Carlo data may by estimated by augmenting with a function that con- 
strains parameters, which are poorly determined statistically, to plausible values. Consider 
for example the logarithm of a Wilson loop, denoted by W ^ which has the expansion 

W{ol) = cia + 020^ + c?,o? + (53) 

We will use ctj [i = 1, 2, . . . , n^) to denote the set of couplings at which the simulations are 
done. We may define an augmented as in Eq. (p!T|), 

X^(Cn) - xL,{Cn) = X\C^) + E (54) 

where the second term on the right tends to constrain poorly determined parameters to the 
range ± o"n, based on our prior experience with the power series. In the examples in this 
Appendix, we will use c„ = and o"„ = 1. 

A well designed calculation should minimize the errors in the final results, for a given 
amount of computer time. The uncertainties in the fit parameters c„ are determined from 
the inverse of the Hessian matrix, which we denote by i^mn; where 

_ 1 

2dCmdCn ^^^^ 



Then the uncertainty in c„ is 



5cn = (56) 
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In the case of Eq. (|53|), an explicit expression for the Hessian matrix can be obtained 

where (J-w{a) is the statistical error in W{a). We note that a useful approximation for the 
statistical error in the Wilson loop (or its logarithm) is 

(Jw{a) ^ /a, (58) 



where / oc 1/V CPU time and is independent of a. 

The optimal selection of the values at which the simulations are to be done may 
be determined by numerically minimizing the 5c„ with respect to the a^. The optimal 
placement of a's for the first couple of parameters may be guessed without doing much 
calculation. For example, since relative statistical errors are independent of a, C\ is obtained 
most accurately by running at the smallest possible coupling a\ (avoiding round-off errors), 
thereby minimizing truncation errors from higher-order terms in the series expansion. The 
smallest total error in C2 is then obtained by choosing a second coupling such that the 
statistical error in the slope of W ja, which is ~ f 012/0.2 = /, is equal to the truncation error, 
which is given by a3al/a2 = (y^ai- Explicit numerical minimization reproduces expectations 
for these simple cases, but also produces optimized design parameters for arbitrary numbers 
of points and allocations of CPU time among them. The results of a three-point optimization 



are illustrated in Fig. |14[ One sees that as the CPU time is increased, smaller and smaller 
couplings are obtained from the minimization calculations, though the optimal a's do not 
fall quite as quickly as the fourth root of CPU time (as would be expected in a two-point 
fit). 

Perturbative series in lattice QCD are typically used in nonperturbative Monte Carlo 
calculations with a's in the range of about 0.15-0.30. To illustrate the effects of optimization 
calculations on the final result of a nonperturbative simulation, we consider the application of 
the perturbative series Eq. ( |53D to a simulation with a = 0.25. The results are shown in Fig. 
|15|. The errors coming from the lowest-order coefficients are greatest when the short-distance 
simulations are done with low statistics, but these errors decrease most rapidly with CPU 
time; hence the order of the term that contributes the largest error to the nonperturbative 
quantity rises as a function of the CPU time of the perturbative simulations. 

Similar results, different in detail, are obtained for other quantities. For example, for the 
third order coefficient K3 for the 5x5 Wilson loop, shown in Fig. ^, the statistical errors are 
//a^ and the truncation error is o"c«. The minimization formula gives the optimal placement 
for a single-point simulation as 

a^{lf, (59) 

as expected. 

Constrained curve fitting formulas for parameter fitting provide a concrete way of trans- 
lating estimates of truncation errors into designs of efficient computations. For one or two 
parameters, they lead to exactly the same the same choices of parameters that intuitive 
guesswork provides. However, they also provide clear optimizations in the much more com- 
mon situations in which we have too many parameters to guess about, or in which we are 
trying to design new runs to improve the results from imperfectly designed initial runs. 
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FIG. 1. Monte Carlo results for for the 5x5 Wilson loop, after the effects of zero modes are 
removed at leading order from the simulation data using Eq. (^). The statistical errors are smaller 
than the plotting symbols. The filled square shows ci from perturbation theory. The dashed line 
shows the results of a fit to Eq. 
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FIG. 2. Monte Carlo results for k^^^ for the 5x5 Wilson loop, after the effects of zero modes are 
removed at leading order from the simulation data. The filled square shows C2 from perturbation 
theory. The dashed line shows the results of a fit to Eq. (^), where ci is constrained to its 
perturbative value. 
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FIG. 3. Monte Carlo results for K2 for the 5x5 Wilson loop, when the effects of the first-order 
zero mode term are not removed from the simulation data. The dashed line shows the results of a 
fit to the data, described in the text. 
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FIG. 4. Results for for the 5x5 Wilson loop. The dashed line shows the results of a fit to 
Eq. (^), where ci and C2 are constrained to their perturbative values. 
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FIG. 5. Simulation results for the temporal Polyakov line on a 4^ lattice at /3 = 9 with periodic 
boundary conditions. Run-time histories are shown for (a) Rc Pt and (b) ImP^. A scatter plot of 
Im Pt versus Re Pt is shown in (c) . 
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FIG. 6. Simulation results for the temporal Polyakov line on a 4^ lattice at /? = 9 with twisted 
Txy boundary conditions. The panels are the same as in Fig. ^. 
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FIG. 7. Simulation results for the temporal Polyakov line on a 4^ lattice at /? = 9 with twisted 
Txyz boundary conditions. The panels are the same as in Fig. ^ 
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FIG. 8. First-order coefficient for tlie tadpole-improved self-energy from perturbation theory, 
using different boundary conditions. The dashed lines show fits to Eq. (^5|). The filled square 
shows the infinite- volume value ci = 1.0701. 
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FIG. 9. Second-order coefficient for the self-energy from perturbation theory, using different 
boundary conditions. The dashed hnes show fits to Eq. (p7|). The filled square shows the infi- 
nite-volume value C2 = 0.117. 
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FIG. 10. Monte Carlo results for ki for the self-energy, Eq. (^91). The results for each lattice 
size L are plotted versus the renormalized coupling ap{q* = 0.84/a). The lowest set of data points 
is for L = 3, and the highest set is for L = 11. The dashed lines show the results of fits to Eq. 
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FIG. 11. First-order coefficient for tfie self-energy from Monte Carlo simulations {<^^) and 
analytic perturbation theory {0^"^). The filled square in (a) shows the perturbation theory value 
of ci on an infinite lattice. The difference between the Monte Carlo results and the perturbation 
theory is shown in (b). 
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FIG. 12. Second-order coefficient for the self-energy from Monte Carlo simulations (c^^*^) and 
analytic perturbation theory (C2"'"). The filled square in (a) shows the perturbation theory value 
of C2 on an infinite lattice. The difference between the Monte Carlo results and the perturbation 
theory is shown in (b). 
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FIG. 13. Simulation results for the third-order coefficient for the static-quark self-energy, after 
subtracting logarithms at 0(1/L), according to Eq. (|50|). The dashed line shows the result of 
a fit to Eq. (|5l|), with the shaded area corresponding to the 68% confidence level region for the 
infinite- volume coefficient. 
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FIG. 14. Values of the couplings 02 (lower curve) and (upper curve) which minimize uncer- 
tainties in the fit parameters ci, C2 and C3, as functions of the CPU time (in arbitrary units). The 
coupling ai is always about zero. 
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FIG. 15. Contribution of errors in the perturbative coefficients to the final error in a 
long-distance simulation of a Wilson loop, as a function of the short-distance simulation time 
(in arbitrary units) that determined the c„. The contribution of a given coefficient c„ to the final 
error is 5c„a", where here we take a = 0.25. The errors coming from the lowest-order coefficients 
are greatest for low statistics in the short-distance simulations, but decrease most rapidly with 
CPU time. 
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